function B = free_slip_bc(V,T,TE,ET,dofs,d,cr,desc,bdr)
% function B_out = outflow_bc(V,T,TE,ET,dofs,d,cr,desc,bdr_out)
% this function implement the outflow boundary condition
d_max = max(d);nbdr = length(bdr); rows = nbdr*(2*d_max-1);
Indx1 = zeros(6*rows,1); % there are 3 nonzeros each row
Indx2 = zeros(6*rows,1);
S = zeros(6*rows,1);
row_begin = 1;
pos = 1;
% figure;hold on;
% begin to treat boundary condition one by one bound edge
for k = 1:nbdr
    eg = bdr(k);  tri = ET(eg,1); dk = d(tri);  
    v1_loc = find(TE(tri,:)==eg);
    v2_loc = mod(v1_loc,3)+1;
    v3_loc = mod(v2_loc,3)+1;
    V1 = V(T(tri,v1_loc),:); 
    V2 = V(T(tri,v2_loc),:);  
    V3 = V(T(tri,v3_loc),:);  
    tri_dofs = get_tri_dof(dofs,tri);
    if v1_loc==2;
        v1_loc = -2;
    end
    % get the dofs on first 3 layers adjacent to boundary edge, which was
    % assumed to be the edge opposite to the first logical node
    line1 = cr_indices(0,dk-1,v1_loc,cr);
    line2 = cr_indices(0,dk-2,v1_loc,cr);
    tau = (V3-V2)/norm(V3-V2);
    % in tcord, the vertices order must use physical order
    A = [1,1,1;V(T(tri,1),1),V(T(tri,2),1),V(T(tri,3),1);...
        V(T(tri,1),2),V(T(tri,2),2),V(T(tri,3),2)];
    t = A\[0;tau(1);tau(2)];
    mat_t = desc_mat(dk,t(1),t(2),t(3),desc);
%     normal = [tau(2),-tau(1)];
    n = A\[0;tau(2);-tau(1)];
    mat_n = desc_mat(dk,n(1),n(2),n(3),desc);
    mat_nn = desc_mat(dk-1,n(1),n(2),n(3),desc);
%     should use this, but the right hand side is 0,so use the convienence
%     Dn = d(tri)*mat_n(line1,:);
%     Dnt = d(tri)*(d(tri)-1)*mat_nn(line2,:)*mat_t;
    Dt = mat_t(line1,:);
%     Dt = mat_t(line1,:);
    Dnn = mat_nn(line2,:)*mat_n;
    [i,j,s] = find([Dt;Dnn]);
    L = length(i);
    Indx1(pos:pos+L-1) = row_begin + i -1;
    tri_dofs = get_tri_dof(dofs,tri);
    Indx2(pos:pos+L-1) = tri_dofs(j);
    S(pos:pos+L-1) = s;
    row_begin = row_begin + 2*dk - 1;
    pos = pos + L;
end
% hold off; axis equal; plot_t(V,T(1:posT,:));
N = row_begin - 1; pos = pos - 1;
dim = max(max(dofs)) - min(min(dofs)) + 1;
B = sparse(Indx1(1:pos),Indx2(1:pos),S(1:pos),N,dim);

function desc = desc_mat(d,lam1,lam2,lam3,desc_pattern)
m_rows = d*(d+1)/2; m_cols = (d+1)*(d+2)/2;
I = (1:m_rows)';       Id = ones(m_rows,1);
desc = sparse(I,desc_pattern(I,1),lam1*Id,m_rows,m_cols) + ...
        sparse(I,desc_pattern(I,2),lam2*Id,m_rows,m_cols) + ...
        sparse(I,desc_pattern(I,3),lam3*Id,m_rows,m_cols);